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The efficient generation of meshes is an important component in the numerical solution of 
problems in physics and engineering. Of interest are situations where global mesh qual¬ 
ity and a tight coupling to the solution of the physical partial differential equation (PDE) 
is important. We consider parabolic PDE mesh generation and present a method for the 
construction of adaptive meshes in two spatial dimensions using stochastic domain decompo¬ 
sition that is suitable for an implementation in a multi- or many-core environment. Methods 
for mesh generation on periodic domains are also provided. The mesh generator is coupled 
to a time dependent physical PDE and the system is evolved using an alternating solu¬ 
tion procedure. The method uses the stochastic representation of the exact solution of a 
parabolic linear mesh generator to find the location of an adaptive mesh along the (artificial) 
subdomain interfaces. The deterministic evaluation of the mesh over each subdomain can 
then be obtained completely independently using the probabilistically computed solutions 
as boundary conditions. The parallel performance of this general stochastic domain decom¬ 
position approach has previously been shown. We demonstrate the approach numerically 
for the mesh generation context and compare the mesh obtained with the corresponding 
single domain mesh using a representative mesh quality measure. 


1 Introduction 

The numerical solution of many partial differential equations (PDEs) benefits from the construc¬ 
tion of an adaptive grid automatically tuned by the solution itself. The quasi-Lagrangian (QL), 
r-refinement, approach used here keeps the number of mesh points and the mesh topology fixed, 
moving the mesh continuously in time using a moving mesh PDE (MMPDE). The solution of 
the mesh PDE gives a continuous mesh transformation between an underlying computational 
co-ordinate and the required physical co-ordinate. Both the mesh and the solution are obtained 
at each time step. The QL approach can be implemented in either an alternating or simulta¬ 
neous manner. The simultaneous QL approach treats the MMPDE and physical PDE as one 
large coupled system. At each time the new mesh and new solution on that mesh are found 
concurrently. Hence the mesh reacts instantly to changes in the physical solution. This highly 
nonlinear coupling may destroy exploitable structure which exists in the discretization of the 
physical PDE alone. The alternating approach uses the current mesh and physical solution to 
update the mesh alone, this new mesh then facilitates the computation of the updated physical 
solution. This introduces a time lag in the mesh as the new mesh is based only on the current 
physical solution. Computationally, however, this decoupling reduces the size of the discrete 
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problem. Furthermore, the solver becomes more modular; the mesh and physical solvers can be 
called in alternating fashion; each solver designed to take advantage of the structure inherent 
in each subsystem. The simultaneous approach is generally thought to be more difficult and 
expensive to solve and hence the alternating method (or a variant thereof) is typically used in 
two or more spatial dimensions. As we will see below, the alternating approach fits well with 
the stochastic domain decomposition approach we describe to parallelize our computations. 

The general QL approach has shown great promise in recent years, solving problems in 
meteorology [7], relativistic magnetohydrodynamics [14], combustion and convection in a porous 
medium [9], groundwater flow and transport of nonaqueous phase liquids [16], Stefan problems 
[4], semiconductor devices [30], and viscoelastic flows [31], phase change problems [3], multiphase 
flows [26], and low speed viscous flow [18], to name just a few. A thorough overview of PDE 
based moving mesh methods may be found in [15]. 

Recently, motivated by the alternating solution method, one of the authors has studied the 
parallel solution of the nonlinear MMPDE alone using a Schwarz based domain decomposition 
approach. In [12], Haynes and Gander propose and analyze classical, optimal and optimized 
Schwarz methods in one spatial dimension at the continuous level. A numerical study of clas¬ 
sical and optimized Schwarz domain decomposition for 2D nonlinear mesh generation has been 
presented in [13]. In [10], a monolithic domain decomposition method, simultaneously solving 
a linear mesh generator coupled to the physical PDE, was presented for a shape optimization 
problem. The authors used an overlapping domain decomposition approach to solve the coupled 
problem. 

In this paper, we present an efficient, parallel strategy for the solution of the moving mesh 
PDE based on a stochastic domain decomposition method proposed by Acebron et al. [1]. The 
motivation is two-fold. First, we wish to reduce (by parallelization) the potential burden of 
having to solve an additional (mesh) PDE. Second, it is often mesh quality, not an extremely 
accurate solution of the mesh PDE, which is important. As we will see, the stochastic domain 
decomposition approach is a means to these ends. 

A stochastic domain decomposition (SDD) method to find adaptive meshes for steady state 
elliptic problems was presented in [5]. The SDD approach, as originally formulated in [1], uses 
a probabilistic form of the point-wise solution for linear elliptic boundary value problems. The 
point-wise solution is evaluated only at the introduced subdomain interfaces. The approximation 
of the solution at each interface point is obtained independently using Monte-Carlo simulations 
and these evaluations are then used as Dirichlet boundary conditions for the (deterministic) 
subdomain solves which can be computed in parallel. The mesh PDEs are generally not solved 
to high accuracy. Mesh quality, allowing an accurate representation of the physical PDE, is 
what is generally required. The parallel algorithm and lower accuracy requirement makes the 
proposed SDD method computationally attractive. The lower accuracy requirement allows one 
to terminate the Monte-Carlo simulations well before convergence. 

The existence of a stochastic representation of the exact solution of linear parabolic problems 
allowed us to extend the SDD approach to (linear) parabolic mesh generators in [5]. There we 
considered the time relaxed form of the Winslow-Crowley variable diffusion mesh generation 
method, first described in [29]. Only the solution of the mesh generator for a specified analytic 
mesh density function was considered. 

Here we consider the generation of time dependent meshes where the mesh PDE is coupled 
to a physical PDE of interest. The coupli ng to the physi cal solution u , is provided by, for 
example, an arc-length type function p = wl + ct(v%. + u*). Furthermore, due to a stochastic 


2 



representation for the solution of PDEs subject to periodic boundary conditions we give, for 
the first time, a method to generate meshes for periodic problems using a stochastic domain 
decomposition approach. 

In Section 2 we describe the time dependent mesh generator used in this paper and how the 
coupling between the physical PDE and mesh PDE is handled for the global, single domain 
solution. Furthermore, we describe a mesh generator for a spatially periodic problem. Section 3 
provides some background on the stochastic solution of linear time dependent, periodic and non¬ 
periodic PDEs and describes how this is used to generate a non-iterative domain decomposition 
algorithm. We precisely illustrate how to obtain approximations to the point-wise solution of 
the mesh PDEs using the stochastic approach and finally how to generate the meshes using the 
stochastic domain decomposition framework. In Section 4 we illustrate the algorithm for various 
examples including Burgers’ equation in both the non-periodic and periodic situations and the 
shallow water equations on a periodic domain. We conclude with some observations and items 
for further study in Section 5. 


2 Mesh generation approach 

As discussed extensively in the references above, in ID the QL r-refinement approach generates 
a physical mesh x E [0,1], via a continuous mesh transformation x(£), where £ is an underlying 
computational variable. A guiding principle is the equidistribution principle of DeBoor (in 
ID) [8,11,28], which finds a mesh transformation x(£) by enforcing 


f£ f 1 

/ p(t,u,x)dx = £ / p(t,U)X)dx, 

J o J o 


where the mesh density function p gives a measure of the level of difficulty or error in the 
physical solution u. The physical solution u may be given analytically or as the solution of a 
physical PDE. In differential form, the mesh transformation may be found as the solution of the 
quasilinear BVP 


d ( , X dx\ 


rr(0) = a, x(l) — b. 


a) 


Here we have suppressed the t and u dependency in p. This gives the physical co-ordinates 
Xi — x(^i) as a function of a (typically uniform) computational grid £^. This BVP can be 
written as a linear BVP in terms of the inverse mesh transformation £(x) as 


d 

dx 


1 d£ 

p(x) dx 


= 0, £(a) = 0, £(&) = 1. 


( 2 ) 


The linearity of the BVP for £(x) makes it easier to solve (in some sense) and in higher 
dimensions has the additional benefit that it is easy to say concrete things about the well- 
posedness of the mesh transformation. As will be discussed below, the linearity of a mesh 
generator is also a prerequisite for the stochastic domain decomposition method. The obvious 
disadvantage is that the solution of the BVP for £(x) does not give the mesh locations in 
the physical variables x. One alternative is to transform the physical PDE from the physical 
variables to the new computational co-oordinates. Alternatively, inverse linear interpolation 
could be used to find the x locations by projecting £(x) onto a uniform £ grid. 

Indeed, in [5] we constructed our stochastic DD method for steady state mesh generation 
using the natural 2D extension of the linear BVP (2). In the time dependent case considered in 
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this paper we will construct a stochastic DD method directly in the physical co-ordinates. The 
time stepping will provide a natural linearization as we will show below. 

A natural way to derive the BVPs above (and which extends to higher dimensions) is as the 
Euler-Lagrange equations whose solutions minimize certain functionals of the required mesh 
transformations. For example minimizing the functional 

iix] =ll ( pw !) 

leads to the quasi-linear BVP (1) above. A minimizer £(x) of the functional 



satisfies the linear BVP (2) above. 

For time dependent PDEs it is useful to use a formulation which involves the mesh speed, 
such a mesh equation is called a moving mesh PDE (MMPDE). If p — p(t,x) the functional /[£] 
becomes 



As described in [15], the direction £ which reduces /[£] is given by the gradient flow equation 

_P d / 1 d£\ 

dt r dx \p{t,x) dx) 

where r > 0 is a user specified constant which determines the response of the mesh to changes in 
p (or u). Here P is a positive definite differential operator that we choose with some flexibility. 

We can switch the roles of x and £ in the gradient flow equation to get a mathematically 
equivalent moving mesh PDE in the physical variables x: 

dx 1 dx f dx\~ 2 f dx\~ l d ( dx\ 

m = rYp{ p ^J (aej 

Hence this resulting equation would retain the nice well-posedness properties. 

If we choose P = (px^) 2 , then we get 

dx 1 d ( dx\ 

m = rHYM)' () 

which is referred to as MMPDE5. Our focus here is the generation of time dependent meshes 
by using this nonlinear parabolic mesh generator. 

The coupling to the physical PDE and physical solution u is provided by the mesh density 
function p(x,u). In a typical deterministic implementation, the mesh and physical solution 
are updated by discretizing MMPDE5 and the physical PDE in time and either solving the 
resulting large nonlinear system of algebraic equations for both the mesh and physical solution 
simultaneously or proceeding in an alternating fashion. The alternating or MP approach [15] 
used in this paper freezes p in the time discretized mesh equation at the current u n to compute 
the next mesh x n+1 and then integrates the physical PDE, using mesh x n+1 , to obtain u n+1 . 
The stochastic solution representation given below requires the PDE to be linear. This MP 
procedure effectively linearizes the MMPDE. 
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In two spatial dimensions we use the mesh generator 

Xt = + V|x, y t = • V^y + V|y, 

P z P ? (4) 

x(0, £) = x 0 = y(0, £) = y 0 = L y y, 

where £ = (£, 77 ) and = (d/d^d/dy), which we solve over the square computational domain 
Ll c — [0,1] x [0,1]. For the sake of simplicity we assume a rectangular domain Ll p = [0 ,L X \ x 
[0, L y \, where L x > 0, L y > 0. At the actual boundaries of the computational domain fi c , we 
employ the fixed boundary conditions x(0,rj) = 0, x(l,rj) — L X: y(£, 0) = 0 and y(£,1) = L y . 
The remaining boundary conditions, rr(£, 0), #(£, 1), y( 0, rf) and y( 1, rf) are found by solving the 
respective one-dimensional forms of the mesh generator (4) along these boundaries. 

Note that each of x and y in mesh generator (4) may be obtained by first dividing the two 
dimensional version of (1) by p and then relaxing the equation in time. We have found in 
practice that this form of a relaxed mesh generator gives better meshes than by relaxing the 
original version (1). In the discrete formulation, these two possible time relaxed mesh generators 
are related by a scaled time step. 

We are also interested in mesh generation on periodic domains. It has been pointed out 
in [27], that in order to use the mesh generator (4) on a rectangular periodic domain of physical 
dimensions Q p — [ 0 , L x [x [0, Ly [, 

x{t, l,r/) = x(t, 0, Tj) + L x: y(t , £, 1) = y(t , f, 0) + L y , 

one should express the mesh and physical PDEs in terms of the displacements X = x — £, where 
£ = (L x £, LyTj), x = (x,y) and X = (X, Y), which are directly periodic on Lt c . In this new set 
of variables, the mesh generator (4) becomes 

x t = ^ • wx + v?x + Y t = ^ + vIy + 

p ^ p p ^ p 

An alternative to the procedure proposed in [27] works directly in the original variables x and 
y and properly extends the computational domain using the periodicity of the grid. This allows 
one to evaluate the derivatives of x on the boundaries without having to take into account the 
actual size of the physical domain Q p . As this approach allows one to work with the physical 
coordinates x directly, we use this second method in this paper. 

3 Stochastic domain decomposition and mesh generation 

In this section we describe stochastic domain decomposition methods for mesh generation on 
periodic and non-periodic domains and describe an implementation which couples the mesh 
generator to the solution of a physical PDE. 

3.1 Stochastic domain decomposition for linear PDEs 

The main idea of stochastic domain decomposition as proposed in [1] (see also [2]) rests on 
the stochastic representation of the exact solution to linear elliptic (resp. parabolic ) boundary 
value problems. This now classical connection between stochastic analysis and boundary value 
problems was first uncovered by Kakutani [19, 20] for the Dirichlet problem using Brownian 
motion. The monographs [21,24] provide a more recent exposition. 
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Numerically evaluating this stochastic representation of the exact solution of a linear bound¬ 
ary value problem using Monte-Carlo methods enables one to compute the point-wise numerical 
solution to the underlying PDE. From the practical point of view, this is fundamentally different 
to solving a PDE using, say, finite differences, which requires the computation of the numerical 
solution over the entire domain even if it is only needed in a single point. 

Notoriously, Monte-Carlo methods converge slowly, with convergence rates proportional to 
TV -1 / 2 where N is the number of Monte-Carlo simulations, if pseudo-random numbers are 
used [25], and hence they play a role mostly for higher dimensional problems, where they can 
be shown to outperform deterministic methods. 

An alternative is to use them in the context of domain decomposition. Namely, splitting the 
entire domain into non-overlapping subdomains, the stochastic solution can be used to com¬ 
pute the point-wise interface solutions between the subdomains. Once these interface solutions 
are determined with sufficient accuracy, they act as Dirichlet boundary values for the single 
subdomains. The PDE solution over each subdomain is computed deterministically using an 
appropriate discretization of the underlying PDE. The main advantage of this method is that 
iteration, as is required in classical domain decomposition methods (such as Schwarz methods), 
can be completely avoided. Also, the solutions over each subdomain can be obtained in parallel 
and thus the method is suitable for massively parallel computing architectures. 

In [5,6] we have shown that the stochastic domain decomposition technique is an effective 
way for the parallel generation of adaptive meshes. A main motivator for the approach is that 
it is, in general, not necessary to compute the meshes with high accuracy. What is important is 
to obtain meshes with high mesh quality. It was shown in [5,6] that even meshes that are not 
accurate solutions to the mesh PDEs can have good mesh quality. This characteristic of mesh 
generation enables an increase in the efficiency of the stochastic domain decomposition method. 


3.2 Stochastic analysis for Dirichlet boundary value problems 

In this section we give the specifics of the stochastic analysis required to generate the solution 
of linear parabolic boundary value problems. Specifically, we consider system (4) where the 
unknowns x — x(t,^rf) and y — y(£,£, 77 ) are required on the time-space domain given by 
[0 ,T] x D c , with T being some finite final time. System (4) is supplemented with the boundary 
conditions x(t,£,rj) |^q c = an d ?/(£,£, ?7) |ao c — ?7), for given continuous functions 

/ and g. The initial values are x(0, £, rf) — xo(£ 5 v) an d 2/(0, £, rf) — yo (£, rf). 

It is important to stress here that (4) is nonlinear if the mesh density function is a function 
on the physical domain, that is p — p(t,x,y). However, as was indicated in Section 2, in the 
practical implementation it is possible to freeze p at time layer t n when computing the mesh at 
time t n+1 . This effectively boils down to a linearization of the mesh generator (4). It is thus 
appropriate to assume that p in (4) is not a function of x and y but of some auxiliary variables x 
and y (and time), i.e. we assume that p = p(t,x,y). Then this mesh generator becomes linear 
and allows for a stochastic representation of its exact solution [23], given by 


x(t,£.,v) = E 

^0 (t ) ) 1 [ran c >t] 

+ E 

f(t - T d n c ,$(T d n c )) l[ran c <i] 

= e 

yo(®(t))l [Tdslc>t ] 

+ E 

g(t - Tqq c > ^ ( r dQ c )) 1 [r S n c <t] 


where the stochastic process <&(£) satisfies the stochastic differential equation (SDE) 


d $(i) = Iv^di + \/2dW(t). 


(5b) 
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In (5), E[-] is the expected value, tqq c is the first exit time of the stochastic process <&(£) starting 
at (£, 77 ), W is two-dimensional Brownian motion and 1 is the indicator function. See [21,23,24] 
for a more extensive discussion of this subject. We note that the stochastic solution (5a) has 
contributions from both the specific initial condition and boundary values of the mesh generator. 

3.3 Stochastic analysis for periodic problems 

The mesh generator (4) is in the form of a system of linear, second order, parabolic PDEs. In 
Section 3.2 we wrote the stochastic point-wise solution assuming a Dirichlet boundary value 
problem. On periodic domains, the celebrated Kac-Feynman formula can be used to obtain the 
stochastic representation of the mesh generator (4), see e.g. [23]. In this case, only initial values 
rr(0, £, 77 ) = xq (£, 77 ) and 7/(0, £, 77 ) = yo(€,v) are given. The solution to (4) can then be written 
as 


= E[rc 0 (^)], y(t, £, rj) = E [y 0 ($)], (6) 

where <l> satisfies the same stochastic differential equation as given in (5b). 

3.4 Stochastic domain decomposition on periodic and non-periodic domains 

While the stochastic solutions (5) and (6) can in principle be used to obtain the solution to 
the parabolic mesh generator at any time t > 0, our mesh generation problem is slightly more 
complicated. The mesh density function p is linked to the solution of the physical PDE system, 
which changes over time. For this reason, in practice we use (5) and (6) only to advance the 
mesh over a single time step, At, from t n to t n+1 . 

We discretize the SDE (5b) using the Euler-Maruyama method, 

$ fc+1 = + YiE A t s + V / 2A4N(0,1), (7) 

P (*»,**) 

with constant time step A t s = A t/K, k = 0,... ,77 — 1, where N is a two-dimensional vector 
of Gaussian distributed random numbers with zero mean and variance one. In other words, one 
time step of size At is split into K sub-time steps. This splitting is necessary so that the use 
of an excessively small time-step At for the solution of the physical differential equation can be 
avoided. The mesh density function p remains fixed at time step t n . The derivatives V^p are 
approximated with finite differences. 

In practice, the starting points of the stochastic process at time £ n , 4>°, coincide with the 
grid points where the stochastic solution is required. Once the new values at time t k+l are 

computed, both p and V^p have to be approximated at 4> /c+1 . For this, bi-linear interpolation 
is used. The procedure is repeated until <f> K is computed, which coincides with the value of 
the stochastic process at time t n+1 . We then evaluate the initial values of xo and yo given at 
time t n at the new location <fr K using bi-cubic interpolation. This gives the values xo(<f> K ) and 
po($ K ), which approximate the actual values xo(<& n+1 ) and yo(<& n+1 ) that are needed in both 
the solutions (5a) (for the first term) and (6). 

Solving Dirichlet boundary value problems stochastically, a boundary test has to be applied 
to determine whether the process has left the domain within the sub-time step [t k ,t k+1 ]. 

A linear Brownian bridge is used as interpolating process as discussed in [17]. If the process left 
the domain before time £ n+1 , the integration can be stopped and the second term in (5a) can 
be evaluated. No such boundary test is needed for periodic domains. 
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In order to estimate the expected value using the arithmetic mean, the procedure is repeated 
N times. 

3.5 Local subdomain problem 

In the domain decomposition solution to the problem (4), we only use the stochastic solutions (5) 
and (6) to generate the subdomain interface solutions for the Dirichlet and periodic problems. 
Once these solutions are computed, the solution to the mesh generation problem on the individ¬ 
ual subdomains becomes a Dirichlet boundary value problem. In order to solve this problem, the 
original parabolic mesh generator (4) has to be solved with the Dirichlet boundary conditions 

x lc>^ = A 

where dQ l 2 c denotes the boundary of the z-th subdomain of the computational domain and f l 
are the values obtained either from the stochastic representation (5) or the boundary conditions 
on the global mesh where d£l l c p| <9D C ^ 0. 

For the local subdomain solver we discretize (4) using centered finite differences for the spatial 
derivatives and a forward Euler method for the time stepping. 

3.6 Domain decomposition solution 

The domain decomposition strategy relies on combining the stochastic evaluation along the 
interfaces with the deterministic subdomain solves. At each time step, the stochastic solution 
procedure provides the Dirichlet boundary conditions for the deterministic single-domain solver. 
This allows the computation of the new mesh at the time step t n+1 . The solution values can 
either be evaluated stochastically at each (<^, rjj) along the artificial interfaces or a sample can 
be evaluated with the rest obtained by interpolation. An optimal placement strategy for the 
interpolation nodes was presented in [5]. 

Once the new mesh is computed, the mesh density function p(t n+1 , x, y) is evaluated on the 
new mesh and the solution procedure is repeated with the new mesh density function. 

3.7 Mesh quality 

As mentioned previously, there is a crucial difference between the SDD method as proposed in [1] 
for the computation of numerical solutions of general linear elliptic boundary value problems and 
for the mesh generation case. The latter does not necessarily require a very accurate solution of 
the mesh equation; here, mesh quality is more important. 

There are several ways of assessing the quality of an adaptive mesh and different mesh quality 
measures based on properties such as equidistribution and alignment have been derived. See [15] 
for an extensive discussion of mesh quality measures. In [5] it was proposed to use the geometric 
mesh quality measure for the assessment of the quality of meshes generated using the SDD 
method. The geometric mesh quality measure is defined by 

1 tr(J T J) 

2 y / det(J T J) ’ ^ ' 

where J is the Jacobian of the transformation x = #(£, 77 ), y = y(£, 77 ) and K is a mesh element 
in D c . The quantity Q(K) measures how far the mesh cell is from being equilateral, that is we 
have Q(K ) > 1 with Q(K) — 1 precisely for an equilateral cell. 





In [5] we argued that using this measure is appropriate as meshes computed using stochastic 
methods show several kinks when they are far away from convergence. Thus, the values of Q(K ) 
are in general larger for such meshes compared to grids that are computed using deterministic 
methods. Of course, the absolute value of Q(K) depends on the mesh density function used and 
the underlying problem for which an adaptive grid is computed. We are therefore not interested 
in the absolute value of Q(K) but only in the ratio between Q SD (W), the geometric mesh quality 
measure of the reference mesh obtained on a single domain, and Q SDD (K ), the geometric mesh 
quality measure of the grid obtained using the SDD method. If this ratio R(K ) = (5 SD /Q SDD 
is close to one, the mesh obtained from the SDD method is a good approximation to the single 
domain mesh. 

The quantity R is a function of each individual mesh cell. In the next section, for the sake 
of convenience, we only list the maximum and mean values, i? max and i? mean respectively, over 
all mesh cells. 

4 Numerical Results 

In this section, numerical experiments in 2D are presented to demonstrate the effectiveness of 
the domain decomposition algorithm and its suitability for problems with both Dirichlet and 
periodic boundary conditions. 

4.1 Burgers’ equation with Dirichlet boundary conditions 

The first test problem considered is the scalar form of the two-dimensional Burgers’ equation 

Ut + (iL 2 /2) x + (u 2 /2) y + v(u xx + Uyy) = 0, (9) 

[x,y] £ = [o, l] x [0,1]. The initial and Dirichlet boundary conditions are chosen such that 

the exact solution is 

“=( 1 + exp ( £± lr ^)) ■ 

and we consider the case of a moderately small diffusion coefficient v — 0.005. In this problem 
the smaller zq the more convection dominates, and large gradients develop and move to the 
boundaries for t > 0, thus requiring a higher concentration of mesh points to resolve the oblique 
shock that propagates with time. This is a common test problem in the moving mesh literature 
[22,32]. The coordinate transformation (x, y) — (x(£, t), y(£, t)) is used to rewrite the 2D Burgers’ 
equation in QL form in computational coordinates, and this is discretized in the computational 
domain using centered differences in space and a trapezoidal rule for the time integration. This 
is coupled to the mesh generator (4). Both (4) and (9) are then solved alternately in time (the 
MP procedure) for each subdomain using an arc-length mesh density function 

P= y/l + a(ul + ul), 

with a = 10/ \\u 2 + UyWoo. The Dirichlet boundary conditions (at the subdomain interfaces) are 
found by evaluating (5a) using (7). At the physical boundary a ID version of (4) is solved. 
We integrate the finite difference discretization of (9) using 41 x 41 grid points on the physical 
domain D c = [0,1] x [0,1] with a time step At = 0.001. The initial time is to = 0.25. We 
use N = 10000 Monte-Carlo simulations with a time step A t s = At/10 for the solution of 
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(a) t=0.75 


(b) t=1.00 


(c) t=1.25 
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Figure 1. Evolution of mesh for Burgers’ equation using the stochastic DD method with 2x2 (top) and 4x4 
(bottom) subdomains. 


the stochastic differential equation (7). In Fig.l the evolution of the mesh is shown when the 
physical domain is split into 2x2 and 4x4 subdomains. 

We report the mesh qualities and the Zoo-norm comparing the exact solution to the numerical 
solution at times tf = 0.75, tf = 1 and tf = 1.25 in Table 1. The geometric mesh quality 
measure is computed for the single domain solution and compared to the DD solution obtained 
by splitting the physical domain into 2 x 2, 3 x 3 and 4x4 subdomains. 


Table 1. Mesh quality and /oo-errors for Burgers’ equation with Dirichlet boundary conditions. 


tf 

/ SD 

6 oo 

/2x2 

L oo 

d 2x2 
^max 

d 2x2 
■^mean 

/3x3 

6 oo 

td 3x3 

1 bn ax 

d 3x3 
Jl mean 

/4x4 

oo 

d4x4 

^max 

d4x4 

■^mean 

0.75 

0.027 

0.023 

0.99 

0.99 

0.031 

0.95 

0.98 

0.026 

0.95 

0.98 

1 

0.031 

0.033 

0.99 

0.99 

0.320 

0.97 

0.99 

0.034 

0.97 

0.99 

1.25 

0.030 

0.031 

1 

l 

0.35 

0.92 

1 

0.24 

l 

1 


Table 1 shows that all the Zoo-errors on meshes generated using SDD are approximately equal 
to the errors found on the associated single domain meshes. Also, the ratios of the geometric 
mesh quality measures of the single domain and DD solutions are close to one for all times, 
indicating that the DD solutions are a good approximation to the single domain solution. 

4.2 Periodic Mesh Generation 

In order to demonstrate the generation of meshes over periodic domains we first study the case 
of a prescribed mesh density function. In particular, we re-visit the five ring problem in the form 
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considered in [5]. That is, we consider an analytically specified velocity function of the form 

2 / \ 2 


u = tanh 

+ tanh 

+ tanh 


R ( x 2 + y 2 — i 


+ tanh 
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+ tanh 


R \ ( x + l) + 


on the domain D p = [—1 , l[x [—1,1[, with R = 30. A simple arc-length mesh density function 


P= + «*), 

with <a = 0.2 is used. The physical domain is discretized using 41 x 41 grid points, the final 
integration time is t — 0.05 with a time step At = 5 x 10 -4 . Since this time step is quite small, we 
could use A t s = A t for the solution of the discretized SDE (7) as well. Again, N = 10000 Monte- 
Carlo simulations were used. In Fig. 2 the mesh generated on a single domain is compared to 
the mesh generated by splitting the physical domain into 4x4 subdomains. 




(a) single domain solution 


(b) solution with 4x4 subdomains 


Figure 2. A comparison of the single domain solution and that obtained for 4x4 subdomains for the five ring 
problem. 


On inspection these meshes appear to be of similar quality and the computed measures of 
mesh quality confirm this, with i? max = 1 and i? mean = 1 for the 2x2, 3x3 and 4x4 
configurations. The mesh quality is identical to the single domain reference case. 


4.3 Burgers’ equation on a periodic domain 

We repeat here the integration of Burgers’ equation (9) but now using a doubly periodic domain 
of size Op = [0, 27 t[x [0, 2n[. The initial condition is 

u = uo + A sin(x + y — 2 ti) , 

where uq = 0.75 and A = 0.5. The diffusion coefficient is v — 0.001. Again, 41 x 41 grid points 
are used and the time step of the integration was At = 0.005, which coincides with the time 
step used in the SDE (7). Only N = 1000 Monte-Carlo simulations were used. The results of 
the integration at times t = 0.85 and t — 1 are collected in Table 2. 

The results in Table 2 confirm that the domain decomposition solution leads to meshes with 
almost the same geometric mesh quality as the single domain reference solution. 
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Table 2. Mesh qualities for Burgers’ equation with periodic boundary conditions. 


tf 

d 2x2 
^max 

d2x2 

^mean 

d3x3 

^max 

d3x3 

■^mean 

d4x4 

■^max 

d4x4 

■^mean 

0.85 

0.99 

l 

l 

l 

0.99 

l 

1 

0.99 

l 

0.99 

l 

0.99 

l 


4.4 Shallow Water Equations on a Periodic Domain 

In this final example we solve the system of shallow-water equations in nondimensional form 
ut + uu x + vu y + h x = 0, 

Vt + uv x + Wy + hy = 0, (10) 

h t + uh x + vh y + h(u x + v y ) = 0, 

where (iq v ) is the two-dimensional velocity field and h is the height of a water column over a 
constant reference level. For this problem, periodic boundary conditions are considered. 

We discretize system (10) in computational coordinates using centered differences for the 
spatial derivatives and a trapezoidal rule for the time integration. The physical domain is of 
size = [0, 27t[x [0, 27t[, which is discretized using 41 x 41 grid points. The time step of the 
integration was At = 0.005, which again coincides with the time step used for the solution of 
the SDE (7). We found that using N — 1000 Monte-Carlo simulations gives sufficiently good 
results. 

The initial condition is a pile of water of height h — 12.5 in the center of the domain over 
a base level at height h = 10. This initial condition simulates the breaking of a dam, i.e. the 
evolution of the water level once the walls holding the pile of water have been removed. 

Since the initial pile of water decays rapidly into gravity waves, we kept the integration times 
short and report our mesh quality results only at t = 0.05 and t = 0.15. A comparison of the 
mesh for the single domain solution and that computed for 4 x 4 subdomains can be seen in 
Fig. 3. 




(a) single domain solution 


(b) solution with 4x4 subdomains 


Figure 3. A comparison of the single domain solution and that obtained for 4x4 subdomains for the shallow 
water equations at t = 0.15. 


As we saw in the previous examples the meshes appear to be of similar quality and this is 
substantiated by the computed mesh quality measures which can be found in Table 3. The mesh 
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qualities obtained from the domain decomposition solution are very close to that of the single 
domain reference solution. 


Table 3. Mesh quality for the shallow-water equations with periodic boundary conditions. 


tf 

d 2x2 
^max 

d 2x2 
^mean 

d 3x3 
Jl max 

d 3x3 
Jl mean 

d4x4 

^max 

d4x4 

^mean 

0.05 

0.99 

l 

0.97 

l 

0.99 

l 

0.15 

0.99 

l 

0.99 

l 

0.99 

l 


5 Conclusion 

Originally, the SDD method has been proposed in [1] for the parallel solution of linear elliptic 
boundary value problems. In [5, 6] we have extended the method to the parallel generation of 
adaptive meshes for prescribed mesh density functions. In this paper we have for the first time 
demonstrated the use of stochastic domain decomposition for the generation of adaptive moving 
meshes for time dependent PDE problems. Furthermore, we have provided the first stochastic 
domain decomposition method for the generation of meshes on periodic domains. 

Our algorithm hinges on the alternating MP procedure for PDE based mesh generation in 
higher dimensions. A fully parallel algorithm would result if the physical PDE was also solved 
in parallel using domain decomposition or other approaches. 

Indeed the linearization provided by the MP procedure has the undesirable affect of decou¬ 
pling the physical solution from the mesh. This effect can be reduced using a M k P procedure. 
This produces a mesh which more closely satisfies the equidistribution principle. In this case, we 
approximate x n+1 by a sequence of sub-meshes x n+1,k , where x n+1,k+1 is obtained from x n+1,k 
by using a step of A t n /K with a linearized MMPDE. In this case p is approximated by p n+1,/c 
obtained by constructing a piecewise linear interpolant of (x n ,p n ) values onto x n+1,k . A M V P 
algorithm, which updates from t n to t n +1 using variable time steps can also be used. The time 
lag between the mesh and physical solution can be reduced by iterating between the mesh and 
physical PDE l times, resulting in a (MP) 1 algorithm. The algorithm ( MP)°° is equivalent to 
the simultaneous solution (if the iteration converges). These solution variants are discussed at 
length in [15]. The implementation of these variants and a study of their efficacy within the 
stochastic domain decomposition framework is a topic of current investigation. 

During our experiments we also detected that fewer Monte-Carlo simulations are needed to 
generate quality periodic meshes. This seems to be due to the absence of an exit time test for 
the periodic case. Work to fully understand and utilize this is also underway. 
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